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' Abstract A Monte-Carlo approach to solving a stochastic jump transition model 

for active-region energy (Wheatland and Glukhov, Astrophys. J. 494, 1998; 
Wheatland, Astrophys. J. 679, 2008) is described. The new method numerically 
solves the stochastic differential equation describing the model, rather than the 
' equivalent master equation. This has the advantages of allowing more efficient 

, numerical solution, the modelling of time-dependent situations, and investigation 

of details of event statistics. The Monte-Carlo approach is illustrated by appli- 
cation to a Gaussian test case, and to the class of flare-like models presented in 
IWheatlandl (|2008p , which are steady-state models with constant rates of energy 
supply, and power-law distributed jump transition rates. These models have 
two free parameters: an index (S), which defines the dependence of the jump 
I transition rates on active-region energy, and a non-dimensional ratio (r) of total 

flaring rate to rate of energy supply. For r <C 1 the non-dimensional mean energy 
(E) of the active-region satisfies (E) ^ 1, resulting in a power-law distribution 
, of flare events over many decades in energy. The Monte-Carlo method is used to 

' explore the behavior of the waiting-time distributions for the flare-like models. 

CNJ , The models with S ^ are found to have waiting times which depart signifi- 

\l ' cantly from simple Poisson behavior when {E) 3> 1. The original model from 

, I Wheatland and Glukhovl l|1998p . with 5 = (no dependence of transition rates 

CNJ ' on active-region energy), is identified as being most consistent with observed 

! flare statistics. 

0\ . 
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5^ ■ 1. Introduction 



Solar flares are explosive events in the solar corona, involving the release of en- 
ergy stored in active-region magnetic flelds. Active regions are dynamic, evolving 
in time due to the emergence and submergence of magnetic flux from the sub- 
photosphere, stressing by photospheric motions, and the occurrence of flares. It 
is of interest to understand the dynamical energy balance of active regions. 
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Flares exhibit a wide range of energies. The largest flares may involve the 
release of up to 10^^ J of energy, and are associated with large-scale expulsions 
of material from the corona (Coronal Mass Ejections, or CMEs), involving com- 
parable energies. Small flare-like events are observed down to the limits of obser- 
vation and are difficult to distinguish from a variety of small-scale solar activity. 



The distribution of flare energies follows a power-law distribution ( Hudson, 1991 



Crosby, Aschwanden, and Dennis, 1993 Aschwanden, Dennis, and Benz, 1998). 
Speciflcally, the frequency-energy distribution J\f{E), i.e. the number of flares 
observed per unit time and per unit energy {E), obeys 



M{E) = AE-^, 



(1) 



where the factor A is a (time-dependent) measure of the total flaring rate, and 
7 w 1.5. This distribution is often constructed for flares from all active regions 
present on the Sun over some period of time, but it also appears to apply to 



individual active regions ( [Wheatland, 2000| ), which suggests that the power law 
is fundamental to the flare mechanism. A popular model explaining the power 



law is the avalanche model (Lu and Hamilton, 1991 Charbonneau et al., 2001), 



in which the magnetic fleld in the corona is assumed to be in a self-organized 
critical state, and subject to avalanches of small-scale reconnection events. The 
distribution ([1]) must have an upper roll-over, to ensure that the total energy 
released in flares is finite ( [Kucera et ai, 1997| [Hudson, 2007P . 

The mechanisms causing flares are not well understood, and flares appear to 
occur randomly in time, although certain properties of active regions correlate 
with flaring (e.ff. [Sammis, Tang, and Zirin[[MinilMcIntos"lil[TMl[Georgoulis and Rust[ | 
120071 [Schrijver[ [^007p . Correspondingly, the prediction of flares is in its infancy: 
the methods used are probabilistic and are rather inaccurate at predicting the 
occurrence of large flares, which are rare but which strongly influence our local 
space weather ( [Wheatland, 2005] [Barnes et al, 2007| [Barnes and Leka, 2008P . 

The occurrence of flares in time may be investigated via a second observable 
distribution, the flare waiting-time distribution, or the distribution of times 
between events (this distribution is also referred to as the "interval distribu- 
tion"). Determinations of flare waiting-time distributions have given varied re- 



sults ( [Pearce, Rowe, and Yeung, 1993[[B"iesecker, 1994[ [Wheatland, Sturrock, and McTiernan, 1998H 



Boffetta et ai, 1 999; Lcprcti rCarbone, and Veltri, 2001] [Moon et al, 2001[ [Wheatland, 2001H 



Moon et ai, 2002j |Kubo, 2008[ ) , suggesting that the observed distribution may 



depend on the particular active region, on time, and that it also may be influ- 
enced by event definition and selection procedures ( Wheatland, 200 1[ Buchlin, Galtier, and Velli, 2005H 
[Paczuski, Boettcher, and Baiesi, 20"05l[Baiesi, Paczuski, and Stella, 2006 ). For some| 
active regions, the distribution is consistent with a simple Poisson process, 
i.e. independent events occurring at a constant mean rate ( [Moon et ai, 200"T| ), 
and the corresponding waiting-time distribution is exponential. Other active 
regions exhibit time-variation in the flaring process, and flare occurrence may 
be approximated by a piecewise-constant, or more generally time- varying Pois- 
son process, in which case the distribution is a sum or integral over expo- 
nentials ( [Wheatland, 200l| ). On longer time scales, the distribution exhibits 
a power-law tail for events from the whole Sun ( [Boffetta et ai, 1999[ ). Some 
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authors have also argued that the process is fundamentally non-Poissonian {e.g. 
[Lepreti, Carbone, and Veltri[ I200T|) . although the arguments neglect the role of 
time-dependence. 

The energy balance of active regions presents a puzzle, because of the large 
drops in energy due to large solar flares. In the following we consider a simple 
"black box" approach to modelling the free energy of an active region (the energy 
available to power flares), an approach that goes back to IRosner and Vaianal 
P978p . Energy is assumed to be continuously supplied to the active region 
by some external mechanism. The energy is stored locally in the corona, and 
some is released at particular times in flares. The flaring process is considered 
to be stochastic, whereas energy supply is deterministic. Flares are treated as 
point processes in time, i.e. they occur at one instant in time, and they involve 
jump transitions (discontinuous changes) in energy. A general model of this 
kind was presented in IWheatland and Glukhovl (|1998p and further developed in 
IWheatlandl (|2008p . The general model is based on a master equation for the 
probability distribution [P{E, t)] for the free energy {E) of an active region at 
time t. The free parameters in the model are a prescribed energy supply rate 
to the system [/3(i?,t)], and prescribed transition rates [a(i?, £", t)], describing 
the rate of jumps from energy E to energy E' < E. These two rate functions 
may have any functional form. In lWheatland and Glukhovl (jlQQSp it was argued 
that the energy-supply rate should not depend on the energy (E) of the system, 
since active regions are driven externally, and hence a constant energy supply 
rate is appropriate. It was also argued that, to produce an appropriate power- 
law flare frequency-energy distribution, transition rates of the form a{E, E') ~ 
{E ~ E')~'' are required (in the steady state). IWheatland and Glukhovl (|1998|) 
and IWheatlandl pOOSp investigated the model by solving the master equation 
in the steady state, for these "flare-like" choices for the free parameters. In 
IWheatland and Glukhovl p998p the emphasis was on the basic model and the 
arguments for the appropriate flare-like choices. A model was constructed that 
could reproduce power-law behavior in the flare frequency-energy distribution 
over an arbitrary number of decades in energy, up to a high energy roll-over set 
by the decline of P{E) at large energy. Hence the flare-like model was confirmed 
to reproduce this aspect of flare statistics. In lWheatlandl (|2008|) the waiting-time 
distributions for the model were also considered, using theory for jump transition 
processes presented for the first time by [Daly and Porporato| ([20071 . It was found 
that the IWheatland and Glukhovl (|1998p model produces an essentially Poisson 
(exponential) waiting-time distribution. A modified model was also considered, 
involving an a{E,E') ~ E{E — E')~^ form for the transition rates. This model 
also reproduced the power-law frequency-energy distribution, but exhibited some 
departure from a simple Poisson waiting-time distribution. 

Master equations may be represented by an equivalent stochastic differential 
equation ( [van Kampen, 19"92{ [Gardiner, 2004P , which provides a complementary 
approach to the problem at hand. Stochastic DEs are amenable to solution by 
Monte-Carlo methods, which in general are simple to numerically implement. 
This paper describes a Monte-Carlo approach to solving the stochastic model 
for solar active-region energy presented in IWheatland and Glukhovl (|1998p and 
IWheatlandl (|2008p . The Monte-Carlo approach has specific advantages over the 
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master-equation approach: it is computationally more efficient, it permits more 
general modelling, in particular the modelling solution of time-dependent prob- 
lems, and it generates an ensemble of flare events and hence permits detailed 
investigation of event statistics. 

The layout of the paper is as follows. The master-equation approach of lWheatland and Glukhovj 
(|1998p and lWheatlandl (|2008p is briefly reiterated in SectionHTTl and the flare-Hke 
choices for the model are explained in Section [^T^ The stochastic DE approach 
to the problem is then presented in Section 13.11 and illustrated by application 
to a Gaussian test case fSection l3.2p . and to the flarc-like cases from IWheatlandl 
pOOSp (Section Is. 3p . including comparison of a Montc-Carlo solution with direct 
numerical solution of the master equation. Section presents a Monte-Carlo- 
based investigation of the variation of the waiting-time distribution for the 
flare-like models, and Section [4] presents conclusions. 



2. Master Equation Approach 
2.1. GENERAL APPROACH 

To begin we briefly reiterate the master-equation formulation of the model, 
following IWheatland and Glukhovl (fTOMP and IWheatlandl (|2008l) . The energy 
[E ~ E{t)] of an active region is assumed to be a stochastic variable which 
evolves in time due to deterministic energy input at a rate f3{E,t), as well as 
due to jumps downwards in energy (flares) at random times and of random sizes, 
described by transition rates a{E,E\t). These are the rate for jumps per unit 
energy from E to E' at time t. The probability distribution [P{E, t)] for the 
energy of the system is given by the solution to the master equation 

= -^[P{E,t)P{E,t)]- \{E,t)P{E,t) 

+ / P{E\t)a{E',E,t)dE\ (2) 
Je 

where 

\{E,t) = / a{E,E',t)dE' (3) 
Jo 

is the total rate of flaring at time t, assuming the system has energy E. [A time 
dependence has been included in the transition rate, by contrast with lWheatland] 
(|2008p .] Two other quantities of interest are the mean total rate of transitions 
(in the average over energy) 



(A) = / X{E,t)PiE,t)dE (4) 
Jo 

and the mean energy of the system 

(E) = / EP{E,t)dE. (5) 
Jo 
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As noted in Section[l] two observable flare distributions are the flare frequency-| 
energy distribution and the waiting-time distribution. The model frequency- 
energy distribution is given by 

poo 

JV{E,t)= P{E',t)a{E',E' - E,t)dE'. (6) 
Je 



|Daly and Porporato (|2007p showed how to obtain the distribution of waiting 



times (r) for jump transition models in the steady state {d/dt = 0). In lWheatlancH 
(|2008p that theory was applied to Equation ^ to yield the model waiting-time 
distribution Pt{t). The details of the derivation are given in IWheatlandl (|2008p . 

In IWhcatla nd and Glukhovl ([T9981) and IWheatlandl (j2008l) the master equa- 
tion was numerically solved in the steady state, for flare-like choices of P{E) 
and a{E,E') (the choices are explained in Section . The methods of solu- 
tion involved discretising the energy as a set of values Ej {i = 1, 2, 3, iV), 
in which case the master equation represents a system of N linear equations 
in N unknowns Pi = P{Ei). Solution of the linear system was performed 
either by relaxation ( [Wheatland and Glukhov, 1998D , or by back substitution 



( Wheatland, 2008 ) . One disadvantage of these methods is that the energy may 



span many decades in energy, in which case a large value of N is required. 
2.2. FLARE-LIKE CHOICES 



In lWheatland and Glukhovl (|1998p . the master equation was solved in the steady 
state {d/dt = 0) for the choices (3{E) = /3o, a constant, and 

a{E, E') = ao{E - E'yOiE - E' - Ec), (7) 

where ao is a constant, Ec is a low-energy cutofl', and 9{x) is the step function. 
In IWheatlandl (pOOg]) Equation ^ was generalised to include an additional 
dependence on the initial energy E: 

a{E, E') = aoE\E - E')-^e{E - E' - E^), (8) 

where (5 is a constant. The problem was solved for Equation ([8]) with (5 = and 
(5=1, and for I3{E) = Pq. 

The physical motivations for these choices is briefly mentioned in Section [1] 
but is worth discussing in more detail. Concerning the energy-supply rate, the 
physical aspect is that the rate does not depend on the energy of the system. This 
is appropriate for a system that is externally driven, and for which there is no 
back reaction of the system on the driver. The picture for the Sun is that energy 
supply comes from below {i.e. from the sub-photosphere), via photosphcric flows 
which cause new fields to emerge into the corona, and twist existing coronal 
fields. The rate at which this occurs is determined by fiow patterns in the sub- 
photosphere. The sub-photosphere is very dense, so it is unlikely that the corona 
can influence the rate of supply of energy, assuming this picture of energy supply 
is correct. Of course, the energy supply rate may depend on time, and the choice 
of a constant supply rate is a simplification. We will return to the question of 
time-dependent driving in Section |4l 
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Concerning the functional forms for the transition rates, first note that sub- 
stituting Equation ([8]) into Equation ([6|) leads to the flare frequency-energy 
distribution 

roc 

Af{E)=aoE-^ {E'f P{E')dE' , (9) 

J E 

for E > Ec-li follows that the frequency-energy distribution is a power law with 
index 7 up to energies E at which P{E) becomes very small. This is consistent 
with the observed power-law frequency-energy distribution ([T]), and the physical 
requirement that the frequency-energy distribution rolls over at large energies 
(to ensure the total mean rate of energy release in flares is finite). An estimate 
of the energy for departure from power-law behavior is provided by the mean 
energy, which may be approximated by ( [Wheatland and Glukhov, 1998D 

/ 2 - 7 \ 1/(^+2-7) 

■ (10) 

In principle, other functional choices leading to power-law behavior are possible, 
although it has proven difficult to identify alternative solutions with an energy- 
supply rate independent of energy that produce a power-law form for J\f{E). 

Another motivation for the choices ([7]) and ([5]) for the transition rates comes 
from consideration of avalanche type models (Lu and Hamilton, 1991 Charbonneau et ai, 200l] ).[ 



In these models the volume involved in flaring is the set of unstable sites which 
trigger one another during the flare "avalanche." The volume of this region 
is found to be scale free, i.e. power-law distributed, and is fractal in shape. 
Assuming the volume of the region is proportional to the energy released, this 
implies a form ~ {E — E')~'^ for the probability per unit time of a transition from 
energy E to E' , assuming flares occur at a constant rate per unit time. Hence 
the flare-likc choices for the master equation may correspond to the avalanche 
model, although the detailed relationship between the two pictures remains to 
be worked out. 

Some support for these choices is provided by the resulting waiting-time 
distributions. The numerical solutions for P{E) in IWheatlandl pOOSp lead to 
waiting-time distributions which are approximately exponential. This may be 
understood by noting that substituting Equation ([8]) into Equation ([3]) gives the 
total rate for flaring 

A(E,<)^{°»^'(^-'"-^"^">*^->'l^^^^" (11) 

For E :$> Ec wc have 

XIE) « -^E~-<+^E\ (12) 
7- 1 

For 6 — 0, Equation implies X{E) is constant (independent of E). For 5 = 1 
we have X{E) oc E, and hence the mean rate may be approximately constant, 
provided P{E) is non-zero only over a fairly limited range in E. Hence the total 
rate of flaring may be approximately constant for 6 = and 6 = 1, consistent 
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with a simple Poisson process. As discussed in Section [U this is compatible with 
observed flare statistics in active regions for which the rate does not vary in time 
dMoon et ai, 2001[ [Wheatland, 200l| . 



3. Stochastic DE Approach 
3.1. GENERAL APPROACH 

Following |Daly and Porporato] (|2007p , the master Equation ^ is equivalent to 
the stochastic differential equation 

^ = f3{E,t)-A{E,t) (13) 

where 

N{t) 

A{E,t) = Y,^E,S{t~U) (14) 

describes the loss in energy due to flaring, with 5{x) being the delta function, 
N{t) being the number of events up to time t, and with the event times ti defined 
by a "state-dependent" Poisson process with occurrence rate X{E, t). The jump 
amplitudes Ai? follow the distribution h{AE,E,t), defined by 

a{E, E - AE, t) = \{E, t)h{AE, E, t), (15) 

so that 

[ h{AE,E,t)d{AE) = 1. (16) 
Jo 

The ODE (fT3|) may be solved in the following way. First, choose a start-energy 
Eg at time t^. The energy of the system evolves deterministically from this time 
up until the first jump at time te = ts + r, where t is a waiting time. The waiting 
time corresponds to a time-dependent Poisson process with a rate A [E{t), t\. To 
evaluate this rate, note that the energy during the deterministic trajectory obeys 

dE 

— =P[E,t). (17) 

Solving Equation ([TT]) with the initial condition E = Eg aX t = tg defines the 
time history E*{t) for the energy, and this together with A = A(i?, t) defines the 
rate A [E*{t), t] prior to a jump. A waiting time may be generated for this rate 



by finding the root of the monotonic function (Wheatland and Craig, 2006): 



fts+T 

F{t) ^\n{l-u)+ I \[E*{t),t]dt, (18) 
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where u is a uniform deviate (a uniformly-distributed number in the range < 
It < 1). If numerical root finding is required, then Newton- Raphson is a suitable 
method (e.g. Press et al. 1992), and in that case it is worth noting that 

F'{t)^\[E*{U),Q, (19) 

where = ts + t. 

Once the waiting time r has been generated, the jump may be simulated. 
The energy before the jump is Ef. = E*{te), the end-energy of the deterministic 
trajectory. The size AE of the jump may be determined by generating a random 
variable from the distribution h{AE, Ee,te), which is defined by Equation (fT5| . 
A value AE may be obtained by the usual technique of transforming a uniform 
deviate to a random variable from the required distribution ( [Press et ai, 19921 ). 
Once AE' is calculated, a new start energy Eg = E^ — AE is specified at time tg 
just after the jump, and the whole process may then be repeated. This procedure 
can be repeated an arbitrary number of times, to give a simulated time history 
of energy E(t) for the system over an arbitrary number of jumps. Provided 
Equations ([T7)) and ^TE\\ are straightforward to evaluate, the process involves 
relatively little computational expense. 



3.2. GAUSSIAN TEST CASE 



To illustrate the method, consider the "Gaussian" test case discussed in lWheatland and Glukhovj 
([T9981) aud i Wheatla^ ([2008]) . namely the case with I3{E, t) = (3o and a{E, E' , t) =| 
ao (where ao and /3o are constants) . In that case the solution to the steady state 
master equation is 

P{E) = aEe-^''^\ (20) 

with a = ao/Po, and from Equation ([6|) the frequency-energy distribution for 
jumps is also a Gaussian: 

MiE)^aae~i''^\ (21) 

From Equation ([3]) the total rate of events is X{E) = aoE, and from Equation ^ 
the mean total rate is (A) = (ao/^o)^^^- Using the method outlined in lWheatlandl 
(|2008p . the waiting-time distribution is also a Gaussian: 

p^^^^^^2a^Y\-^.„0or^^ (22) 

To simulate the Gaussian test case using the Monte-Carlo approach, note that 
the solution to Equation (fT7|) with /3{E, t) ~ /?□ and with starting energy Es at 
time tg is 

E*{t)^Eg+po{t~tg), (23) 
and we have X{E) = a^E, so 

X[E*{t)] = ao[Eg+f3o{t-ts)]. (24) 



SOLA: MCActiveRegioiiEnergyArXiv.tex; 3 February 2009; 2:39; p. 8 



Monte-Carlo Simulation of Active Region Energy 



9 



Equation p8|) evaluates to 



F{t) = ln(l - u) + aoEsT + -Q^o/3oT^ 



(25) 



and taking the positive root of this quadratic function gives 




1 + 



aE^ 1 - u 



) 



1/2 



1 



(26) 



For this case the distribution of jump amplitudes, from Equation (jlSp . is given 



for < AE < E. The jump energies are uniformly distributed on (0, E), and a 
jump may be generated from a uniform deviate u using AE = Eu. Finally, the 
mean energy 



provides a suitable starting energy for simulation. 

Figure [T] illustrates a Monte-Carlo solution of the Gaussian test case using 
Equations ((23l) - (p8|) . for the choices ao = 0.01 and (3q = 1. The upper panel 
in the figure shows the time history of system energy for the simulation for 
25 jumps, and the lower panel shows the corresponding event energies versus 
time. The energy of the system grows linearly with time between jumps, and the 
jumps occur with a rate which increases linearly with energy. The jump sizes 
arc uniformly distributed, up to the current energy of the system. 

Figure [2] illustrates the Monte-Carlo solution of the Gaussian test case with 
the same parameters (ao = 0.01 and /3o ~ 1) for 5000 waiting times and jumps, 
and compares the results with the analytic expressions. The left-hand panel 
shows the histogram of the system energy, and the right-hand panel shows the 
histogram of waiting times. The corresponding analytic distributions P{E) and 
Pt{t), given by Equations ([20]) and ([22]) respectively, are shown by the solid 
curves. The histogram of system energy was obtained by sampling the simulated 
time history of energy E{t) at 5000 random times, uniformly distributed over 
the duration of the simulation. These results illustrate the application of the 
Monte-Carlo approach, and confirm that it reproduces the analytic results. 

3.3. FLARE-LIKE CASES 

Next we consider the flare- like cases corresponding to Equation ([8]) , from lWheatland[ 
pOOSp . For those case the total rate of flaring is given by Equation pT|) . and the 
energy supply rate is /9 = /So, a constant. If the last event is at t = is, when the 
energy is Eg, then the subsequent deterministic trajectory in energy (prior to 
the next jump) is given by Equation Using this result and Equation pT|) . 



by 



h{AE,E) = ^ 



(27) 




1/2 



(28) 
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Figure 1. Monte-Carlo solution of the Gaussian test case with qq = 0-01 a-^d /3o = 1, with 
25 jump transitions and waiting times. The upper panel shows the system energy versus time, 
and the lower panel shows the event energies versus time. 




10 20 30 
Energy 



10 20 30 
Waiting time 



50 



Figure 2. Monte-Carlo solution of the Gaussian test case, involving 5000 jump transitions, 
and parameter values oq = 0.01, /3o = 1. The histogram in the left-hand panel shows the 
distribution of system energy, and the histogram in the right-hand panel shows the distribution 
of waiting times. The corresponding analytic results are shown in the two panels by solid curves. 



Equation ^TE\\ evaluates to 



F{t) = ln(l -u) + 



ao/f3o r E-^+^ r 



7-1 
1 



5-7 + 2 



(29) 
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where e — Es ii Es > Ec, and e ~ Ec ii Es < Ec. In this case numerical root 
finding is required, and for the application of Newton-Raphson it is helpful to 
note from Equation p9p that 



F'(r) = ^ {E, + /3ot)' k-^+i - {Es + /3ot)^^+'' 
7-1 L 

The distribution of jump energies defined by Equation (jlSp is 

(7-l)(Ai?)-^+i 



hiAE, E) 



(30) 



(31) 



for Ec < AE < E. A variable with this distribution may be generated from a 
uniform deviate u for a given E using the transformation 



AE = [E-^+^ - u {E-''+^ - E-'<+^)] 



-7 + l\] -1/(7-1) 



(32) 



A suitable starting energy for a simulation is provided by Equation ([TO]). 

Following IWheatland and Glukhovl ([TOM)) and I Wheatland) ([^005)) . it is useful 
to non-dimcnsionalize, by introducing 



E' 



and 



_ OLaE[ 
r = 



'(5-7+2 



(33) 



(34) 



Equations ()29p and ([50)) become 



F{t) = ln(l -u) + 



7-1 1,(5+1 



{Es+t) 



1 



5-7 + 2 

where 1 = Eg Eg >\, and e = 1 if iJ^ < 1 , and 



F'ir) 



7-1 



(E.^+r)' 



1 - ^B , 



-7 + 1' 



respectively. The transformation used to generate jump energies is 



AE = 

and the starting energy is 



l-u[l-E-'^') 



-1/(7-1) 



2-7 



1/(5+2-7) 



(35) 



(36) 



(37) 



(38) 
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The parameters 6 and r define the specific model being solved. Equation pS]) 
implies that r < ^ is required for (E) > 1. Smaller values of f lead to larger values 
of (E) and hence more decades of power-law behavior in the frequency-energy 
distribution, as explained in Section 12.21 Many decades of power-law behavior 
are observed for flares on the Sun, implying a small value of r. 

As an example of applying Equations ((35|) - (|38p . we consider the case 6 = 0, 
from IWheatland and Glukhovl p?98p . Figure [3] shows the results of a Monte- 
Carlo solution with i5 = and r = 0.02. The upper panel shows the time history 
of the energy of the model active region over the first 50 jumps, and the lower 
panel shows the corresponding flare energies versus time (with a logarithmic 
scale). The mean energy of the system is (E) = 625, which is the starting energy 
for the simulation. Figure [3] illustrates the character of the flare-like models. 
The active-region energy grows linearly with time between flares, flare sizes are 
power-law distributed, and the total flaring rate is approximately constant, since 
E ^ 1. During this period of time only relatively small flares occurred, with the 
largest event having energy close to 100 in non-dimensional units. 




Figure 3. Monte-Carlo solution of the flare-like case from Wheatland and Glukhov (1998), 
with 5 = and r = 0.02. The upper panel shows the active-region energy versus time, and 
the lower panel shows the flare energy versus time, for a period of time including 50 jump 
transitions. 

Figure |4] illustrates the Monte-Carlo solution with the same parameters ((5 = 
and r = 0.02) for 3 x lO"* waiting times and jump transitions, and compares the 
results with direct solution of the master equation. The upper-left panel shows 
the time history of the energy of the system, with the mean energy (E) = 625 
shown by a horizontal line, which is also indicated by an arrow near the left-hand 
axis of the panel. The upper-right panel shows the histogram of the energy of 
the system, obtained by sampling the simulated time history of energy at 3 x 10^ 
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random times. The mean energy is shown by a solid vertical line. The solid curve 
is the distribution obtained by solving the master equation in the steady state, 
using the method in l Wheatlandl (|2008p . The lower-left panel shows the histogram 
of waiting times, in a log-linear representation, together with the waiting-time 
distribution obtained from the solution to the master equation (solid curve). 
The lower-right panel shows the flare frequency-energy histogram, together with 
the distribution obtained from the solution to the master equation (solid curve) , 
and the mean energy (solid vertical line). The Monte-Carlo solutions agree with 
the direct solution of the master equation. These results confirm the finding in 
IWheatlandl (|2008p that the waiting-time distribution is essentially exponential 
in this case. 




Figure 4. Montc-Carlo solution of the flare-like case from Wheatland and Glukhov (1998), 
involving 3 X 10'* jump transitions, and the parameter value r = 0.02. The figure shows the ac- 
tive-region energy versus time (upper left), and histograms of: the active-region energy (upper 
right); the waiting-time distribution (lower left); and the flare frequency-energy distribution 
(lower right). The corresponding results obtained by solving the master equation in the steady 
state are shown by solid curves. The active-region mean energy is indicated by an arrow near 
the left-hand axis of the upper-left panel, and is shown by a solid vertical line in the two panels 
on the right. 

The comparison in Figure |4] illustrates the computational advantage of the 
Monte-Carlo method over direct solution of the master equation. The results 
shown by the solid curves in Figure [¥] require the solution of a linear system 
in 5000 unknowns, to obtain sufficient energy resolution to ensure accuracy. By 
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comparison, the Monte-Carlo solution requires the simulation of 3 x 10** waiting 
times and jump transitions, but each of these calculations is simple. As a result 
the Monte-Carlo method is substantially faster than the linear solution. 

3.4. WAITING-TIME DISTRIBUTIONS FOR FLARE-LIKE CASES 

In IWheatlandl (|2008p it was shown using numerical solutions of the master 
equation that the flare-like models with S ~ and 5 — 1 exhibit approxi- 
mately Poisson (exponential) waiting-time distributions, for certain choices of 
the non-dimensional ratio r = aoE^~^^^ / Po of transition rates to energy-supply 
rate [in IWheatlandl (|2008p this ratio was labelled ocq, or just ao, when the bar 
was dropped]. These results are briefly discussed in Section [221 However, the 
results also showed some departure from the simple Poisson distribution. In this 
section we further investigate the behavior of the models and in particular the 
waiting-time distributions, using Monte-Carlo solutions. 

Figure [5l illustrates three solutions for the case (5 = 0. Each solution involves 
3 X 10** waiting times and jump transitions. The upper row in Figure 5 shows 
the flare frequency-energy distributions in a log-log representation, and the 
lower row shows the corresponding waiting-time distributions, in a log-linear 
representation. The left-hand pair of distributions is for r = 0.5, the center pair 
is for r = 0.05, and the right-hand pair is for r = 0.005. The solid vertical 
lines in the frequency-energy distributions show the approximation to the mean 
energy (E) given by Equation (j38p . The three choices of r shown correspond 
to values (E) = 1, {E} = 10^, and (E) = 10^. The solid Hues on the waiting- 
time distributions show the exponential form A,„e^^'"'^, where Am is the overall 
mean rate of events, i.e. the number of events divided by the total time (non- 
dimensionaliscd). The upper row of Figure 5 confirms that the frequency-energy 
distribution is a power law with index 7 below a roll-over set by the largest energy 
the system is likely to attain, which may be roughly approximated by (E). For 
smaller values of r, the system attains larger energies, as flaring is less frequent. 
The lower row of Figure [5] shows that the waiting-time distribution becomes 
exponential as (E) increases (or as r decreases). This may be understood using 
the argument given in Section l^?^ for f ^ 1, the energy E of the system satisfies 
i? ^ 1, in which case the approximation of Equation (|12p applies, and the total 
rate of flaring X{E) is then independent of E. 

Figure [HI illustrates three solutions for the case 6 = 1, again with 3 x 10^ 
waiting times and jump transitions. The format of the flgure is the same as for 
Figure \S\ and the values of r for the three cases are chosen so that the values 
of {E) [using the approximation of Equation (|38p ] are the same, i.e. (E) = 1 
(left), {E) = 10^ (center), and (E) = lO"* (right). The corresponding values of r 
are given above the three rows in the flgure. The upper row of Figure [6l shows 
the expected power-law behavior below an upper roll-over given approximately 
by (E) . The lower row of Figure [6l shows that the waiting-time distribution is 
approximately exponential for (E) = 10^, but for (E) = 10^ there is an excess of 
large waiting times by comparison with the exponential form. Referring again to 
the argument in Section [2?2l for r <C 1 we expect E ^ 1, in which case the total 
flaring rate varies approximately as X{E) oc E. The linear variation in mean 
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Figure 5. The flare frequency-energy distributions (upper row) and flare waiting-time distri- 
butions (lower row) for the flare-like case with 5 = 0. The left-hand pair of distributions is for 
(E) ?s 1 (r = 0.5), the center pair is for (E) Ri 10^ (r = 5 X 10~^), and the right-hand pair is 
for (E) rj 10"* (r = 5 X 10"^). 




Figure 6. The flare frequency-energy distributions (upper row) and flare waiting-time distri- 
butions (lower row) for the flare-like case with 5 = 1. The left-hand pair of distributions is for 
(E) = 1 (r = 0.5), the center pair is for (E) = 10^ (r = 5 X 10^^), and the right-hand pair is 
for (E) = lO"* (r = 5 X lO"'^). 



flaring rate with energy leads to some departure from the simple exponential 
form, with the departure becoming more significant as r decreases. 

Figure [7| illustrates three solutions for the case (5 = 2, again with 3x10^ 
waiting times and jump transitions. The format of the figure is the same as for 
Figures [5] and [6l with the values of r again corresponding to approximate mean 
energies (E) = 1 (left), (E) = 10^ (center), and (E) = 10* (right). The results 
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Figure 7. The flare frequency-energy distributions (upper row) and flare waiting-time distri- 
butions (lower row) for the flare-like case with 5 = 2. The left-hand pair of distributions is for 
(E) = 1 (r = 0.5), the center pair is for (E) = 10^ (r = 5 X 10"'^), and the right-hand pair is 
for (E) = 10* (r = 5 X 10""). 



for the frequency-energy distribution (upper row) are as expected, with power- 
law behavior below an upper roll-over given approximately by (E). The results 
for the waiting-time distribution are similar to the case 5=1, but are more 
pronounced. For (E) = 10^ there is very approximate exponential behavior, with 
some excess of large waiting times. For {E) = 10^ there is substantial departure 
from the exponential model, with a pronounced excess of large waiting times. 

— 2 

Using the argument in Section [2T21 we expect approximate oc E dependence of 
the total flaring rate for E ^ 1. This variation in the mean flaring rate with 
energy leads to the departure from the exponential form. 

The results shown in Figures O-IT] suggest that the (5 = model may be the 
preferred flare-like solution. The flare frequency-energy distribution is observed 
to be a power law over many decades in energy, which implies r ^ 1 . The models 
with 5^0 exhibit signiflcant departure from Poisson waiting-time statistics for 
r ^ 1. However, it is not clear that any such departure is observed for flares 
on the Sun: some active regions that produce very large flares appear to exhibit 



simple Poisson waiting-time statistics (Moon et at, 2001 Wheatland, 2001) 



The variation in total flaring rate with energy of the system for the cases with 
i5 ^ also implies observable consequences. For example, if we pick out "large" 
events from a simulation ensemble with S ^ and plot the waiting time before 
the event versus the waiting time after the event, then we expect that the waiting 
times after the event will tend to be larger than the waiting times before, since 
large events deplete the system energy, and hence reduce the total rate. Figured 
illustrates this effect. The figure is constructed using the simulation shown in 
Figure [7| with S = 2 and (E) = 10*, and the threshold for a large event is chosen 
to be Q.2{E). The solid line is the line of equality of the two waiting times. As 
expected, the waiting times after the large events tend to be larger than the 



SOLA: MCActiveRegionEnergyArXiv.tex; 3 February 2009; 2:39; p. 16 



Monte-Carlo Simulation of Active Region Energy 



17 



waiting times before the event (more points lie below the line of equality than 
above). No such effect is observed when a similar plot is constructed for the 
model with 5 = and {E) = 10^. This effect is not observed for flares on the 



Sun (Wheatland, Sturrock, and McTiernan, 1998 Wheatland, 2001). 
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Figure 8. Plot of the waiting time before an event, versus the waiting time after an event, for 
all events larger than 0.2(E), for the simulation with 5 = 2 and {E) = 10^ shown in Figure[7] 



4. Conclusions 

This paper presents a Monte-Carlo method for solving the stochastic model for 
active region energy presented in lWheatland and Glukhovl(|1998|) and lWheatland [ 
pOOSp . which in particular is suited to solving the flare-like cases in those 
papers. The method numerically solves the stochastic differential equation de- 
scribing the system, rather than the equivalent master equation, and provides 
a computationally-efficient approach to the problem. The method is demon- 
strated on a simple Gaussian test, and is compared with a direct solution of the 
steady-state master equation for a flare- like case from IWheatland and Glukhovl 

(HMHl). 

The method is used to further investigate the class of flare-like models from 
IWheatlandl (|2008p . which feature constant energy-supply rates /3o and flare 
transition rates of the form 



a{E, E') = aoE^{E - E')-''6{E - E' - E^), 



(39) 



where ao is a constant, and E^. is a low-energy cutoff. The index 7 = 1.5 is the 
observed flare frequency-energy power law index, and (5 is a positive constant 
(we consider the cases 5 = 0, (5 = 1, and 5 = 2). The emphasis in the investiga- 
tion is on the waiting-time distributions for these models, and their adherance 
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to/departure from a Poisson (exponential) form, as a function of S and of the 
dimensionless ratio r = aQE^~'''^'^ / The models require small values off to 
produce flares with a frequency-energy distribution exhibiting a power law over 
many decades [a lower bound to the departure from power-law behavior is set by 
the estimate of the mean energy {E) « [(2 — ^) /rf'^'^^^'^~^^]. For the 5 = Q model 
it is found that the waiting-time distribution becomes a close approximation to a 
simple exponential for small f . This may be explained in terms of the total rate of 
flaring becoming constant for large which applies when f <C 1 . For the 6 = 1 
and 5 = 2 models the waiting-time distribution is approximately exponential for 
intermediate values of r, but exhibits an excess of large waiting times for r ^ 1. 
This result may make the (5 = case the preferred flare-like solution, since it 
is unclear that any such departure is observed for flares on the Sun. Also, the 
dependence of the total flaring rate on the active-region energy for the 5=1 
and 5 = 2 models implies variations in flare rate with flare occurrence which do 
not appear to be observed on the Sun. 

It is interesting to reconsider the correspondence between the stochastic model 
and the avalanche model for flares (Lu and Hamilton, 1991 Charbonneau et ai, 2001] ) ,[ 
in light of the new results. The investigation in this paper excludes the possibility 
of a stochastic model with transition rates of the form of Eq. ([39]) with 5^0, 
which simultaneously exhibits a wide range of flare energies and has uncorrelated 
waiting times. This is due to large flares changing the system energy and hence 
total flaring rate significantly. The choice for the transition rates was motivated 
in part by the avalanche model (as outlined in Section [2T2|l . Avalanche models 
exhibit power-law frequency-energy distributions over many decades, and they 
have uncorrelated waiting times ( Wheatland, Sturrock, and McTiernan, 19"98l ). 
However, in avalanche models the largest flares deplete only a small fraction 
of of the avalanche grid "energy" (see e.g. Fig. 3 in lCharbonneau et fflZ.|[200ip . 
The discrepancy between the two pictures might be due to different definitions 
of energy (the avalanche model "energy" may correspond to free energy plus 
background magnetic energy unavailable for flares). Another possibility is that 
the 5 = model provides the closest match to the avalanche picture. A more 
careful analysis of the correspondence between the two models is needed to 
resolve this point. 

The Monte-Carlo approach has several advantages over solution of the master 
equation. First, as discussed in Section [3.31 the Monte-Carlo solution is numeri- 
cally simple to implement and computationally efficient. Second, as highlighted 
by Figure [51 since the Monte-Carlo approach produces an ensemble of flare 
events, it permits detailed investigation of event statistics. A third advantage 
of the Monte-Carlo method is that it permits solution of the model in time- 
dependent situations. In principle solutions may be constructed for arbitrary 
time variation of the rates P{E, t) and a{E, E' , t). The methods of solution of the 
master equation presented in IWheatland and Glukhovl (|1998p and IWheatlandl 
(|2008p apply only in the steady state, and in particular the method of obtaining 
the waiting-time distribution requires a steady state. As discussed in Section [l] 
time variation in observed flaring rates plays a role in determining the observed 
waiting-time distribution, so it is of interest to consider time-dependent models. 
These models will be addressed in future work. 
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It is possible to construct more general stochastic models for active-region 
energy. Active regions may lose energy by mechanisms other than flaring, e.g. 
due to slow ohmic dissipation of electric currents. If the energy loss is considered 
to be deterministic, and smaller in magnitude than the energy input rate, then 
/3(£^, t) in the present model may be interpreted as a net energy-supply rate, 
and the model may be considered to already accommodate energy losses of this 
kind. If the energy loss is assumed to occur via small random decrements, then a 
suitable model may be to include a Fokker-Planck (diffusion) term in the master 
equation, describing continual small random increases and decreases in energy. 
(The increases may correspond to fluctuations in the energy supply rate.) A 
more general form for the jump-transition master equation, sometimes called 
the Chapman-Kolmogorov equation ( [Gardiner, 2004[ ) includes drift (represent- 
ing deterministic energy input or loss), diffusion (small random input or loss), 
and jump transitions. An active-region model of this form was briefly discussed 
in lWheatlandl (|2008p . In the context of the Monte-Carlo approach, the diffusion 
term corresponds to a Wiener process, and a different method of solution of the 
stochastic DE is then required. 

Models of the kind presented in this paper are difficult to test in detail against 
flare observations because of the difficulties associated with determining the free 
energy of active regions, and the rate of energy supply to active regions. How- 
ever, we note that improved solar observations (and analysis techniques) may 
eventually provide such information. Independent of this, the models provide 
important qualitative checks on our understanding of energy balance in solar 
active regions. For example, as noted in Section 12.21 it has proven difficult to 
identify choices other than the flare-like ones investigated here that produce 
suitable power-law flare frequency-energy distributions. The models are also of 
intrinsic interest because of their description of a dynamical balance involving 
scale-free transitions, and it is possible that they provide suitable descriptions 
of a variety of other physical systems exhibiting power-law behavior. 
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